options(future.globals.maxSize = 891289600)
### Make sure all the packages are installed
if (!requireNamespace("Seurat", quietly = TRUE))
install.packages("Seurat")
if (!requireNamespace("dplyr", quietly = TRUE))
install.packages("dplyr")
if (!requireNamespace("devtools", quietly = TRUE))
install.packages("devtools")
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if (!requireNamespace("RcppML", quietly = TRUE)) {
# BiocManager::install("fgsea", update = FALSE)
# BiocManager::install("limma", update = FALSE)
# devtools::install_github("zdebruine/singlet", upgrade = FALSE)
devtools::install_github("zdebruine/RcppML")
devtools::install_github("zdebruine/RcppSparse")
}
if (!requireNamespace("sparseMatrixStats", quietly = TRUE)) {
install.packages("sparseMatrixStats")
}
if (!requireNamespace("inflection", quietly = TRUE)) {
install.packages("inflection")
}
if (!requireNamespace("org.Hs.eg.db", quietly = TRUE)) {
BiocManager::install("org.Hs.eg.db")
}
if (!requireNamespace("clusterProfiler", quietly = TRUE)) {
BiocManager::install("clusterProfiler")
}
### Load all the necessary libraries
library(Seurat)
library(dplyr)
library(RcppML)
library(RcppSparse)
library(ggplot2)
library(sparseMatrixStats)
library(inflection)
library(glue)
library(tidyr)
library(org.Hs.eg.db)
library(clusterProfiler)
set.seed(687)7 - Factor analysis
Introduction
In this notebooks we are going to carry out Factor Analysis analysis using RcppML. You can see the GitHub repository here. This implementation of NMF is extremely fast and enables us to use large dataset since it works well with sparse matrices. The author is currently implementing it within the Singlet package to work nicely with single cell and soon it will be available!
How does RcppML’s NMF algorithm work preprint
What are the differences between NMF and PCA - This stats.StackExchange summarises quite well what the differences are between NMF and PCA. It uses the image below to visually represent the differences - see the free book An Introduction to Statistical Learning for more details!
Here are some papers that nicely present applicabilities of NMF:
To identify transcriptional programs in our data
To identify shared transcriptional programs across tumors:
- [Hallmarks of transcriptional intratumour heterogeneity across a thousand tumours](10.1038/s41586-023-06130-4)
- [Cancer cell states recur across tumor types and form specific interactions with the tumor microenvironment](10.1038/s41588-022-01141-9)
Key Takeaways
NMF identifies sets of genes “metagenes” representing the main characteristics of the data.
Choosing the rank (K) of the matrix is an important step since it will determine how many “metagenes” are present in our dataset. We can use cross validation to find the best K.
Whe can visualize how important each genes is for each factor (aka - metagene) and how important is each factor for each cell. This way we can determine which metagenes are explaining each cell’s transcriptome.
When looking at metagenes learned across multiple patients and conditions it is imperative to check that the signal is not being provided just by one replicate of our experiment. Not checking this could lead to incorrect interpretation of the data!
Library
Load data
We’re going to be working with a dataset from the paper - Immunophenotyping of COVID-19 and influenza highlights the role of type I interferons in development of severe COVID-19 Download data from cellxgene portal.
# Download the data in data/ directory
# download.file(
# url = "https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds",
# destfile = "../data/workshop-data.rds",
# method = "wget",
# extra = "-r -p --random-wait")
# We can also use the CLI with the wget command below
# wget https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds
se <- readRDS("../data/workshop-data.rds")Set color palettes
disease_pal <- c("#41AE76", "#225EA8", "#E31A1C")
names(disease_pal) <- c("normal", "influenza", "COVID-19")
# flu <- RColorBrewer::brewer.pal(12, name = "YlGnBu")
# normal <- RColorBrewer::brewer.pal(12, name = "BuGn")
# covid <- RColorBrewer::brewer.pal(12, name = "YlOrRd")
donor_pal <- c(
"#66C2A4", "#41AE76", "#238B45", "#006D2C",
"#41B6C4", "#1D91C0", "#225EA8", "#253494", "#081D58",
"#FFEDA0", "#FED976", "#FEB24C", "#FD8D3C",
"#FC4E2A", "#E31A1C", "#BD0026", "#800026")
names(donor_pal) <- c(
"Normal 1", "Normal 2", "Normal 3", "Normal 4",
"Flu 1", "Flu 2", "Flu 3", "Flu 4", "Flu 5",
"nCoV 1", "nCoV 2", "nCoV 3_4", "nCoV 5",
"nCoV 6", "nCoV 7_8", "nCoV 9_10", "nCoV 11"
)Analysis
Convert ENSEMBL IDs to Gene Symbols
Right away we can see how ensembl ids are used in the rownames. Let’s transform them into their matched symbols to make them human-readable:
head(rownames(se))[1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457"
[5] "ENSG00000000460" "ENSG00000000938"
Convert to gene symbols
gene_df <- readr::read_csv(file = "../data/cov_flu_gene_names_table.csv")
symbol_id <- data.frame(index = rownames(se)) %>%
left_join(gene_df) %>%
pull(feature_name)
# re-create seurat object
mtx <- se@assays$RNA@data
rownames(mtx) <- symbol_id
se <- CreateSeuratObject(counts = mtx, meta.data = se@meta.data)Factor Analysis
For the purpose of this vignette we’re going to focus on CD8 and NK cells this object to pull out interferon signalling program:
se <- se[, se$Celltype %in% c("CD8, non-EM-like", "CD8, EM-like", "NK cell")]First we need to normalize the data
se <- NormalizeData(se, verbose = FALSE) %>%
FindVariableFeatures(nfeatures = 3000, verbose = FALSE) %>%
ScaleData(verbose = FALSE) %>%
RunPCA(verbose = FALSE)
VlnPlot(
se,
features = c("CD3D", "CD3E", "CD8B", "NKG7", "FCGR3A"),
group.by = "Celltype",
pt.size = 0,
split.by = "disease") + theme(legend.position = "bottom")Extract the normalized expression matrix and remove genes that are all 0s
library(Matrix)
mtx <- se@assays$RNA@layers$data
colnames(mtx) <- colnames(se)
rownames(mtx) <- rownames(se)
# Remove genes that are all 0s
mtx <- mtx[sparseMatrixStats::rowSums2(mtx) > 0, ]Cross-validation for rank determination
The first step is to determine the rank of our matrix - by this we mean, which is the optimal number of factors we need to decompose our matrix. Here we run cross-validation 3 times across a range of different ranks (k).
start <- Sys.time()
cv_data <- crossValidate(mtx, k = seq(1, 36, 5), reps = 2, verbose = TRUE, tol = 1e-04)
Sys.time() - startTime difference of 4.106416 mins
Now let’s visualize the crossvalidation results
ggplot(cv_data, aes(x = k, y = value, color = rep, group = rep)) +
geom_point() +
geom_line() +
theme_minimal()We can see how it starts to plateau at ~k=20 so we will use that in this script.
Run NMF
k: Number of factors we want to identify
tol: Stands for tolerance and is how small we want the error to be, the smaller the number the better the decomposition but the longer its gonna take to converge. Values in the 10^-5 and 10^-6 return very good results.
L1: Introduces sparsity into the factors and loadings with the aim of removing noisy genes from factors and noisy cells contributing to factors. L1 normalization uses Lasso regularization to increase sparsity and remove the unimportant features.
set.seed(7)
nmf_ls <- RcppML::nmf(
data = mtx,
k = 20,
tol = 1e-06,
L1 = c(0.05, 0.1), # l1 for c(w, h)
verbose = TRUE)Explore the NMF results
W - contains which genes are relevant for each factor. H - contains how important each factor for each cell.
w <- nmf_ls@w
w[1:5, 1:5] nmf1 nmf2 nmf3 nmf4 nmf5
TSPAN6 0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
DPM1 0 9.446519e-05 4.755987e-05 1.160628e-04 7.280245e-05
SCYL3 0 0.000000e+00 0.000000e+00 5.235962e-07 0.000000e+00
C1orf112 0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
FGR 0 5.544984e-05 4.062468e-04 5.537400e-05 2.077896e-04
h <- nmf_ls@h
h[1:5, 1:5] AAACCCAAGAATGTTG-12 AAACCCAAGCTACGTT-6 AAACCCAAGGCCTGCT-12
nmf1 5.462850e-05 5.187584e-05 3.464711e-05
nmf2 7.132754e-05 1.925005e-04 0.000000e+00
nmf3 1.210247e-04 2.170797e-05 0.000000e+00
nmf4 1.459895e-04 0.000000e+00 0.000000e+00
nmf5 5.263035e-05 0.000000e+00 1.276532e-04
AAACCCAAGTGTAGAT-5 AAACCCACACAATGCT-5
nmf1 6.166142e-05 5.885350e-05
nmf2 9.268590e-05 1.126141e-04
nmf3 0.000000e+00 0.000000e+00
nmf4 4.154653e-06 2.351866e-05
nmf5 7.246126e-05 1.722137e-05
Add NMF to Seurat object
se[["FA"]] <- Seurat::CreateDimReducObject(
embeddings = t(nmf_ls$h),
loadings = nmf_ls$w,
assay = "RNA",
key = "Factor_")
se[["FA"]]A dimensional reduction object with key Factor_
Number of dimensions: 20
Number of cells: 19262
Projected dimensional reduction calculated: FALSE
Jackstraw run: FALSE
Computed using assay: RNA
Examine and visualize NMF loadings a few different ways
print(se[["FA"]], dims = 1:5, nfeatures = 5)Factor_ 1
Positive: MALAT1, MT-CO3, MT-CO2, RPL10, EEF1A1
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 2
Positive: GNLY, CD52, GZMH, FGFBP2, NKG7
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 3
Positive: GNLY, GZMB, NKG7, GZMH, IL32
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 4
Positive: GNLY, FGFBP2, NKG7, TYROBP, CCL5
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Factor_ 5
Positive: NFKBIA, GNLY, CD69, CCL4, IER2
Negative: TBCEL-TECTA, RP1-111C20.3, RP5-1067M6.6, H2BE1, RP11-826N14.7
Violin plots for the loadings
VlnPlot(
se,
features = colnames(se[["FA"]]),
ncol = 4,
group.by = "Celltype",
pt.size = 0) + theme(legend.position = "right")VlnPlot(
se,
features = colnames(se[["FA"]]),
ncol = 4,
group.by = "Celltype",
split.by = "donor_id",
pt.size = 0,
cols = donor_pal) + theme(legend.position = "right")We can also visualize this by cell type, condition and patient!
# Preprocess dataset
dd <- bind_cols(se@meta.data, se@reductions$FA@cell.embeddings) %>%
tidyr::pivot_longer(
cols = glue::glue("Factor_{1:ncol(nmf_ls$w)}"),
names_to = "k",
values_to = "loading") %>%
group_by(Celltype, disease, donor_id, k) %>%
summarise(median_loading = median(loading))
lapply(glue::glue("Factor_{1:ncol(nmf_ls$w)}"), function(i) {
dd %>%
filter(k == i) %>%
ggplot(aes(x = disease, y = median_loading, fill = disease)) +
geom_boxplot(outlier.shape = NA, alpha = 0.7) +
geom_jitter(aes(color = donor_id), height = 0) +
facet_wrap(~ Celltype) +
scale_fill_manual(values = disease_pal) +
scale_color_manual(values = donor_pal) +
scale_x_discrete(limits = names(disease_pal)) +
theme_classic() +
labs(title = glue::glue("Factors by condition in {i}")) +
guides(fill="none")
}) %>% patchwork::wrap_plots(ncol = 4)Factors 5 & 7 seem interesting, lets explore them a bit more in depth!
# w[1:5, c("nmf5", "nmf7")]
# Extract top genes from factors 4 & 6 using Unit Invariant Knee
FA_genes <- lapply(c("nmf5", "nmf7"), function(i) {
y_vec <- sort(w[, i], decreasing = TRUE)
# Define inflecion point
n <- uik(x = seq_len(length(y_vec)), y = y_vec)
# Define top 10 genes
df <- data.frame(gene = names(y_vec), value = y_vec, rank = 1:length(y_vec)) %>%
mutate(lab = if_else(rank < 15, gene, NA_character_))
print(ggplot(df, aes(x = rank, y = y_vec)) +
geom_point() +
geom_vline(xintercept = n) +
ggrepel::geom_text_repel(aes(label = lab)) +
labs(title = glue("Factor-{i}")) +
theme_classic())
y_vec[1:n]
})names(FA_genes) <- c("nmf5", "nmf7")
dd <- lapply(names(FA_genes), function(i){
data.frame(
value = FA_genes[[i]],
gene = names(FA_genes[[i]]),
factor = i)
}) %>%
bind_rows() %>%
pivot_wider(names_from = factor, values_from = value, values_fill = 0)
DT::datatable(dd)By looking at these genes we can get a sense of which are the major processes captures by each factor.
NMF-5 seems to be capturing cytotoxic and activation signals since it contains genes encoding for chemokines (CCL4, CCL3), granzymes (GZMB, GNLY, NKG7) and other immune related processes (NFKBIA, CD69…).
NMF-7 contains a wide array of myeloid cell markers, such as S100A8, S100A9, LYZ and genes encoding for MHC-II (CD74, HLA-DRB1, HLA-DRA) as well as cytotoxic genes (GZMA, GZMH, NKG7…). The most straight forward explanation is that these could be doublets or have a high amount of ambient RNA in the soup.
We can take a look at which genes are driving factor 7, Flu-1 donor has very high score compared to the rest and we want to make sure it is not driving that factor by itself. To do so we will select the top 25 genes with the highest loadings in factor 7 and visualize their scaled expression across all donors.
# We'll take a look at the top 25 genes
names(head(FA_genes[["nmf7"]], 25)) [1] "S100A9" "S100A8" "PRF1" "NKG7" "GZMA" "IL32" "CD74"
[8] "GZMB" "LYZ" "GZMH" "HLA-B" "B2M" "HBB" "TXNIP"
[15] "S100A11" "MT-CO1" "PSME2" "PSMB9" "PLAAT4" "RPS4Y1" "CST7"
[22] "TRAC" "FGFBP2" "ACTB" "HLA-A"
# Scale their expression
se <- ScaleData(se, features = names(head(FA_genes[["nmf7"]], 25)))
# Check scale.data
mtx <- se@assays$RNA@layers$scale.data
rownames(mtx) <- names(head(FA_genes[["nmf7"]], 25))
colnames(mtx) <- colnames(se)
# Choose color palette
palette_length <- 100
my_color <- colorRampPalette(c("Darkblue", "white","red"))(palette_length)
# set color breaks
my_breaks <- c(seq(quantile(mtx, .01), 0, length.out=ceiling(palette_length/2) + 1),
seq(0.05, quantile(mtx, .99), length.out=floor(palette_length/2)))
# Note that we are setting the min and max of the scale to .01 and .99 so we exclude extreme values from dampening our signal
# We can see how the max value of the matrix is +10 but q99 is 2
Hmisc::describe(as.vector(mtx))as.vector(mtx)
n missing distinct Info Mean Gmd .05
481550 0 287115 1 -2.226e-06 1.1 -1.59884
.10 .25 .50 .75 .90 .95
-1.31044 -0.59795 0.02677 0.69204 1.16664 1.48214
lowest : -12.5361 -8.98419 -8.74596 -8.11004 -7.98887
highest: 9.87275 9.88233 9.96698 9.97462 10
pheatmap::pheatmap(
mtx,
cluster_rows = TRUE,
cluster_cols = TRUE,
show_rownames = TRUE,
show_colnames = FALSE,
treeheight_col = 0,
annotation_col = se@meta.data[, c("disease", "donor_id", "Celltype")],
annotation_colors = list("disease" = disease_pal, "donor_id" = donor_pal),
color = my_color,
breaks = my_breaks)Carry out GSEA
# read GSEA markers
gsea_ls <- lapply(FA_genes, function(i) {
# http://yulab-smu.top/clusterProfiler-book/chapter5.html#go-gene-set-enrichment-analysis
gsea_results <- clusterProfiler::gseGO(
geneList = i,
ont = "BP",
OrgDb = org.Hs.eg.db,
keyType = "SYMBOL",
minGSSize = 10,
maxGSSize = 300,
pvalueCutoff = 0.1,
pAdjustMethod = "BH",
seed = TRUE)
gsea_results
})
names(gsea_ls) <- names(FA_genes)Visualize top enriched gene sets per cluster
lapply(names(gsea_ls), function(i) {
# Extract gsea
gsea <- gsea_ls[[i]]
gsea <- clusterProfiler::simplify(gsea, cutoff = 0.7)
gsea@result <- gsea@result %>%
dplyr::arrange(dplyr::desc(NES)) %>%
dplyr::filter(p.adjust < 0.1)
tmp_plt <- gsea@result %>%
dplyr::top_n(n = 20, wt = enrichmentScore) %>%
ggplot2::ggplot(.,
ggplot2::aes(
x = NES,
y = forcats::fct_reorder(Description, NES),
size = setSize,
color = p.adjust)) +
ggplot2::geom_point() +
# ggplot2::geom_vline(
# xintercept = 0,
# linetype = "dashed",
# color = "red") +
scale_color_viridis_c(option = "plasma") +
ggplot2::theme_minimal() +
ggplot2::labs(title = i)
}) %>% patchwork::wrap_plots(ncol = 2)And lastly we will save GSEA to a spreadsheet so we can take a look:
gsea_xlsx <- lapply(names(gsea_ls), function(i) {
print(i)
# Extract gsea
gsea <- gsea_ls[[i]]
gsea <- clusterProfiler::simplify(gsea, cutoff = 0.7)
gsea@result <- gsea@result %>%
dplyr::arrange(dplyr::desc(NES)) %>%
dplyr::filter(p.adjust < 0.1)
})[1] "nmf5"
[1] "nmf7"
names(gsea_xlsx) <- names(gsea_ls)
openxlsx::write.xlsx(gsea_xlsx, file = "../data/GSEA.xlsx")Session Info
sessionInfo()R version 4.3.1 (2023-06-16)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] Matrix_1.6-5 clusterProfiler_4.8.2 org.Hs.eg.db_3.17.0
[4] AnnotationDbi_1.62.2 IRanges_2.34.1 S4Vectors_0.38.1
[7] Biobase_2.60.0 BiocGenerics_0.46.0 tidyr_1.3.1
[10] glue_1.7.0 inflection_1.3.6 sparseMatrixStats_1.12.2
[13] MatrixGenerics_1.12.3 matrixStats_1.2.0 ggplot2_3.5.0
[16] RcppSparse_1.0 RcppML_0.5.6 dplyr_1.1.4
[19] Seurat_5.0.2 SeuratObject_5.0.1 sp_2.1-3
loaded via a namespace (and not attached):
[1] fs_1.6.3 spatstat.sparse_3.0-3 bitops_1.0-7
[4] enrichplot_1.20.0 devtools_2.4.5 HDO.db_0.99.1
[7] httr_1.4.7 RColorBrewer_1.1-3 backports_1.4.1
[10] profvis_0.3.8 tools_4.3.1 sctransform_0.4.1
[13] DT_0.29 utf8_1.2.4 R6_2.5.1
[16] lazyeval_0.2.2 uwot_0.1.16 urlchecker_1.0.1
[19] withr_3.0.0 prettyunits_1.1.1 gridExtra_2.3
[22] progressr_0.14.0 cli_3.6.2 spatstat.explore_3.2-6
[25] fastDummies_1.7.3 scatterpie_0.2.1 sass_0.4.8
[28] labeling_0.4.3 spatstat.data_3.0-4 readr_2.1.4
[31] ggridges_0.5.6 pbapply_1.7-2 yulab.utils_0.1.4
[34] foreign_0.8-84 gson_0.1.0 DOSE_3.26.2
[37] parallelly_1.37.0 sessioninfo_1.2.2 rstudioapi_0.15.0
[40] RSQLite_2.3.1 generics_0.1.3 gridGraphics_0.5-1
[43] crosstalk_1.2.1 ica_1.0-3 spatstat.random_3.2-2
[46] vroom_1.6.3 zip_2.3.0 GO.db_3.17.0
[49] ggbeeswarm_0.7.2 fansi_1.0.6 abind_1.4-5
[52] lifecycle_1.0.4 yaml_2.3.8 qvalue_2.32.0
[55] Rtsne_0.17 grid_4.3.1 blob_1.2.4
[58] promises_1.2.1 crayon_1.5.2 miniUI_0.1.1.1
[61] lattice_0.21-8 cowplot_1.1.3 KEGGREST_1.40.0
[64] pillar_1.9.0 knitr_1.45 fgsea_1.26.0
[67] future.apply_1.11.1 codetools_0.2-19 fastmatch_1.1-4
[70] leiden_0.4.3.1 downloader_0.4 ggfun_0.1.4
[73] data.table_1.15.0 remotes_2.4.2.1 vctrs_0.6.5
[76] png_0.1-8 treeio_1.24.3 spam_2.10-0
[79] gtable_0.3.4 cachem_1.0.8 openxlsx_4.2.5.2
[82] xfun_0.42 mime_0.12 tidygraph_1.2.3
[85] survival_3.5-7 pheatmap_1.0.12 ellipsis_0.3.2
[88] fitdistrplus_1.1-11 ROCR_1.0-11 nlme_3.1-163
[91] ggtree_3.8.2 usethis_2.2.2 bit64_4.0.5
[94] RcppAnnoy_0.0.22 GenomeInfoDb_1.36.3 bslib_0.6.1
[97] irlba_2.3.5.1 rpart_4.1.19 vipor_0.4.5
[100] KernSmooth_2.23-22 Hmisc_5.1-0 colorspace_2.1-0
[103] DBI_1.1.3 nnet_7.3-19 ggrastr_1.0.2
[106] tidyselect_1.2.0 processx_3.8.2 bit_4.0.5
[109] compiler_4.3.1 htmlTable_2.4.1 plotly_4.10.4
[112] shadowtext_0.1.3 checkmate_2.2.0 scales_1.3.0
[115] lmtest_0.9-40 callr_3.7.3 stringr_1.5.1
[118] digest_0.6.34 goftest_1.2-3 spatstat.utils_3.0-4
[121] rmarkdown_2.25 XVector_0.40.0 base64enc_0.1-3
[124] htmltools_0.5.7 pkgconfig_2.0.3 fastmap_1.1.1
[127] rlang_1.1.3 htmlwidgets_1.6.4 shiny_1.8.0
[130] jquerylib_0.1.4 farver_2.1.1 zoo_1.8-12
[133] jsonlite_1.8.8 BiocParallel_1.34.2 GOSemSim_2.26.1
[136] RCurl_1.98-1.12 magrittr_2.0.3 Formula_1.2-5
[139] GenomeInfoDbData_1.2.10 ggplotify_0.1.2 dotCall64_1.1-1
[142] patchwork_1.2.0 munsell_0.5.0 Rcpp_1.0.12
[145] ape_5.7-1 viridis_0.6.4 reticulate_1.35.0.9000
[148] stringi_1.8.3 ggraph_2.1.0 zlibbioc_1.46.0
[151] MASS_7.3-60 plyr_1.8.9 pkgbuild_1.4.2
[154] parallel_4.3.1 listenv_0.9.1 ggrepel_0.9.5
[157] forcats_1.0.0 deldir_2.0-2 Biostrings_2.68.1
[160] graphlayouts_1.0.0 splines_4.3.1 tensor_1.5
[163] hms_1.1.3 ps_1.7.5 igraph_2.0.2
[166] spatstat.geom_3.2-8 RcppHNSW_0.6.0 reshape2_1.4.4
[169] pkgload_1.3.2.1 evaluate_0.23 BiocManager_1.30.22
[172] tzdb_0.4.0 tweenr_2.0.2 httpuv_1.6.14
[175] RANN_2.6.1 purrr_1.0.2 polyclip_1.10-6
[178] future_1.33.1 scattermore_1.2 ggforce_0.4.1
[181] xtable_1.8-4 RSpectra_0.16-1 tidytree_0.4.6
[184] later_1.3.2 viridisLite_0.4.2 tibble_3.2.1
[187] aplot_0.2.2 memoise_2.0.1 beeswarm_0.4.0
[190] cluster_2.1.4 globals_0.16.2